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This  study  presents  the  model  development  and  validation  of  an  efficient  method  for  numerical  predic¬ 
tions  of  the  performance  characteristics  of  proton  exchange  membrane  (PEM)  fuel  cells.  To  efficiently 
execute  extensive  modeling  in  a  short  time  period,  the  computational  model  uses  a  half-cell  at  the 
cathode  with  only  one  oxidant  channel  as  the  modeling  domain.  The  original  three-dimensional  (3D) 
model  is  replaced  with  a  reduced  two-dimensional  model.  In  terms  of  these  computational  models  and 
reduced  dimensions,  mathematical  formulations  describing  physical  phenomena  within  fuel  cells  are 
developed.  In  addition,  two-phase  Darcy’s  laws,  which  are  originally  only  applied  to  porous  media,  are 
further  extended  to  describe  the  transport  and  formation  of  liquid  water  within  channels.  In  this  study, 
the  Fortran  programming  language  is  used  as  a  numerical  program  to  implement  iterative  calculations 
and  predict  the  overall  performance  characteristics  for  three  different  cases.  By  comparing  the  modeling 
and  experimental  polarization  curves  for  case  1  and  the  modeling  of  points  of  current  density  vs.  cell 
voltage  ( IV )  of  different  computational  models  for  cases  2  and  3,  this  computational  model  is  found  to 
possess  a  certain  degree  of  accuracy  and  reliability.  In  addition,  the  complete  extensive  modeling  also 
works  more  efficiently  and  with  fewer  computational  resources. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Although  the  use  of  fossil  fuels  drives  economic  development 
and  advancement  of  civilization,  the  environmental  pollution  and 
greenhouse  effects  accompanying  these  fuels  has  encouraged  the 
search  for  cleaner,  more  efficient,  richer  fuels.  Hydrogen  energy 
appears  to  be  a  viable  option.  The  interest  in  hydrogen  energy  has 
recently  led  to  increasing  interest  in  fuel-cell  research  and  applica¬ 
tions.  Proton  exchange  membrane  (PEM)  fuel  cells  are  particularly 
widely  used  with  products  that  have  an  output  power  ranging  from 
1 W  to  1 00  kW  because  of  their  high  efficiency,  low  pollution,  quiet¬ 
ness,  and  low-temperature  start-up  characteristics.  The  future  for 
PEM  fuel  cells  using  hydrogen  energy  is  greatly  anticipated  and  is 
of  interest  to  many  people. 

Research  on  PEM  fuel  cells  involves  fluid  dynamics,  heat  trans¬ 
fer,  electrochemistry,  materials,  structure  analysis,  and  integration 
and  control  of  systems.  Therefore,  harmony  and  coordination 
between  various  technologies  is  needed  to  achieve  stable  and  reli¬ 
able  fuel-cell  performance.  Among  these  technologies,  fuel-cell 
flow  design  plays  a  key  role  such  that  modeling  research  that 
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utilizes  comprehensive  computational  fluid  dynamics  (CFD)  or  sim¬ 
plified  numerical  analysis  has  shown  great  advancements  in  both 
quality  and  quantity.  In  the  1990s,  Springer  et  al.  [1  ]  first  presented 
numerical  modeling  of  PEM  fuel  cells.  Subsequently,  Bernardi  and 
Verbrugge  [2]  and  Nguyen  and  White  [3]  also  developed  fuel¬ 
cell  models.  However,  during  this  period,  limited  computational 
tools  and  resources  resulted  in  modeling  research  focused  on 
one-dimensional  (ID)  models.  In  the  late  1990s,  fuel-cell  model¬ 
ing  gradually  developed  toward  two-dimensional  (2D)  models.  At 
this  stage,  equations  of  continuity,  momentum,  energy,  species, 
current  and  two-phase  flows  were  used  to  describe  the  electro¬ 
chemistry  and  transport  phenomena  within  fuel  cells.  For  example, 
Gurau  et  al.  [4]  evaluated  polarization  curves  under  various  oper¬ 
ating  conditions  with  a  2D  fuel-cell  model  and  analyzed  local 
variations  of  properties,  such  as  oxygen  and  water  vapor  con¬ 
centration  and  current  densities.  Um  et  al.  [5]  investigated  the 
effects  of  hydrogen  dilution  on  cell  performance  and  polarization 
with  a  2D  transient  model.  It  is  worth  mentioning  that  the  addi¬ 
tion  of  two-phase  models  further  improved  the  prediction  ability 
of  fuel-cell  modeling.  He  et  al.  [6]  used  liquid-water  transport 
equations  governed  by  the  shear  force  of  gas  flows  and  capil¬ 
lary  force  to  perform  simulation  and  investigation  for  fuel  cells 
with  interdigitated  flow  fields.  Wang  et  al.  [7]  analyzed  liquid- 
water  transport  phenomena  within  the  cathode  of  fuel  cells  using 
a  two-phase  mixture  model  [8],  in  which  the  effects  of  flood¬ 
ing  on  cell  performance  were  evaluated.  Since  then,  research  on 
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the  2D,  two-phase  modeling  of  fuel  cells  has  been  published 
continually. 

By  2000,  mainstream  research  on  three-dimensional  (3D)  fuel¬ 
cell  modeling,  which  allowed  researchers  to  understand  the 
physical  phenomena  and  mechanisms  within  fuel  cells  more 
clearly,  began  to  develop.  Mazumder  and  Cole  [9,10]  performed 
3D  fuel-cell  simulations  with  gas-phase  and  two-phase  models, 
respectively.  Their  work  not  only  evaluated  the  effects  of  flooding 
on  cell  performance,  but  also  discussed  the  mechanisms  governing 
liquid-water  transport.  Hu  et  al.  [11,12]  built  3D  fuel-cell  mod¬ 
els  with  straight  and  interdigitated  channels,  in  which  the  cell 
performances  of  different  channels  are  compared,  and  transport 
phenomena  and  physical  mechanisms  for  two-phase  flows  are  also 
analyzed  in  detail.  In  addition,  the  model  is  especially  impressive 
for  large-scale  fuel-cell  modeling.  Meng  and  Wang  [13]  built  a  3D 
fuel-cell  model  with  5-channel  serpentine  channels  using  1  million 
grids.  This  was  the  first  time  that  modeling  was  realized  using  the 
parallel  computing  methodology.  Wang  and  Wang  [14]  developed 
a  50  cm2  fuel-cell  model  with  36  channels  and  used  2.7  million 
grids  to  perform  simulations  of  gas-phase  flows.  Then,  Wang  and 
Wang  [15]  also  developed  a  200  cm2  fuel-cell  model,  in  which  only 
gas-phase  flows  were  considered  and  23.5  million  grids  were  used 
in  parallel  computations.  Inoue  et  al.  [16]  performed  simulations 
for  fuel  cells  with  225  cm2  active  areas,  in  which  the  cell  perfor¬ 
mances  of  serpentine  channels  with  different  channel  numbers  and 
patterns  are  compared.  Subsequently,  Inoue  et  al.  [16]  used  their 
previously  developed  model  [17,18]  to  evaluate  the  effects  of  com¬ 
binations  of  different  gas  and  coolant  channels  on  cell  performance. 

After  2007,  the  interest  in  large-scale  fuel-cell  modeling  seems 
to  have  declined.  Instead,  small-scale  fuel-cell  modeling  received 
more  attention.  Most  of  these  studies  evaluated  the  effects  of  flow 
design  on  cell  performance.  For  example,  Wang  et  al.  [19]  investi¬ 
gated  how  straight  and  interdigitated  channels  affect  the  flooding 
and  cell  performances  under  different  operating  voltages.  Kuo 
et  al.  [20]  performed  fuel-cell  simulations  with  straight  and  wave 
channels  to  analyze  microscopic  phenomena  and  compare  overall 
performances.  Min  [21]  presented  fuel-cell  stepped  flow  designs, 
and  evaluated  the  effects  of  the  channel  with  various  step  numbers 
on  cell  performances.  Labato  et  al.  [22]  utilized  CFD  to  predict  the 
performances  of  fuel  cells  with  serpentine,  straight,  and  pin-type 
channels  under  different  operating  conditions,  and  analyzed  how 
these  channel  patterns  affect  the  distributions  of  current  densities. 
In  addition,  a  large  number  of  studies  address  fuel-cell  model 
developments.  For  example,  Sui  et  al.  [23]  built  a  comprehensive 
3D  fuel-cell  model,  and  the  model  was  validated  based  on  the 
comparisons  of  modeling  and  experimental  results  involving 
fuel  cells  with  straight  channels.  Ly  et  al.  [24]  presented  reduced 
fuel-cell  models  to  decrease  the  computational  requirements  and 
preserve  geometrical  resolution,  and  the  modeling  results  are  not 
only  consistent  with  the  experimental  findings,  but  these  models 
also  reduce  the  required  computational  resources  and  time.  Kim 
et  al.  [25]  developed  fuel-cell  models  with  reduced  dimensions 
and  compared  the  modeling  results  of  different  dimensional  fuel 
cells  with  a  single  straight  channel  to  evaluate  which  combinations 
of  dimensions  possess  the  best  efficiency  and  accuracy.  In  addition 
to  the  above  flow  designs  and  model  developments,  there  has 
been  continual  published  research  about  the  modeling  of  fuel-cell 
operating  parameters,  cold  start  and  transient  response.  For 
example,  Yuan  et  al.  [26]  utilized  CFD  to  analyze  the  effects  of 
operating  parameters  such  as  the  outlet  pressure,  cell  temperature, 
relative  humidity  and  stoichiometric  ratio  on  cell  performances. 
Jiao  and  Li  [27]  performed  fuel-cell  simulations  of  various  cold 
start  processes  and  evaluated  the  ways  to  reduce  water  freezing 
during  a  cold  start.  Khajeh-Hosseini-Dalasm  et  al.  [28]  developed  a 
two-phase  transient  model  of  the  fuel-cell  cathode  side  to  investi¬ 
gate  the  time  variation  of  liquid  water  and  gas  transport,  and  good 


agreements  between  modeling  and  experimental  results  verified 
this  model. 

The  decline  in  large-scale  fuel-cell  modeling  research  after  2007 
may  have  occurred  because  this  type  of  work  requires  extensive 
computational  resources,  and  only  a  few  research  institutes  with 
access  to  parallel  computing  can  perform  the  modeling.  It  can  be 
determined  from  the  aforementioned  previous  papers  that  a  sim¬ 
ulation  for  large-scale  fuel  cells  usually  requires  several  hours  to 
several  tens  of  hours,  and  thus,  extensive  resources  and  time  must 
be  devoted  to  conducting  modeling  research.  Although  large-scale 
fuel-cell  modeling  has  received  less  attention  recently  for  these 
reasons,  this  type  of  research  is  often  so  difficult  to  handle  that 
researchers  avoid  it  as  much  as  possible.  The  fuel  cells  used  in  com¬ 
mercial  products  are  usually  large-scale  cells.  Therefore,  in  order 
to  solve  the  above  difficulty,  this  study  presents  a  fuel-cell  model¬ 
ing  method  programmed  in  the  Fortran  computing  language  that 
can  be  applied  to  PEM  fuel  cells  with  various  scales.  This  method 
can  not  only  predict  overall  performance  characteristics,  such  as 
polarization  curves,  reactant  flow  rates  and  pressure  drops,  but  can 
also  provide  local  property  variations  to  further  analyze  physical 
phenomena  and  mechanisms  within  fuel  cells.  The  computational 
models  developed  in  this  study  are  described  in  detail  in  the  fol¬ 
lowing  sections. 

2.  Model  development 

In  this  study,  model  development  covers  three  sections: 
computational  models,  basic  assumptions,  and  mathematical  for¬ 
mulations.  Each  section  is  discussed  separately. 

2.1.  Computational  model 

In  order  to  simplify  fuel-cell  modeling,  a  cathode  physical  model 
that  includes  only  the  regions  covered  by  one  oxidant  channel  is 
considered  in  this  study,  as  shown  in  Fig.  1(a)  and  (b).  In  Fig.  1  (a),  the 
computational  domains  surrounded  by  dotted  lines  in  the  x-y  plane 
include  only  active  areas,  but  they  may  also  include  active  areas  and 
non-active  areas,  both  of  which  are  covered  by  one  oxidant  channel 
and  two  adjacent  half  ribs.  In  Fig.  1(b),  a  cross-sectional  model  in 
the  y-z  plane  is  shown  that  contains  one  oxidant  channel,  adjacent 
ribs,  a  cathode  gas  diffusion  layer  (CGDL),  a  cathode  catalyst  layer 
(CCL),  and  a  half  PEM  for  the  domains,  including  active  areas.  How¬ 
ever,  for  the  domains  of  non-active  areas,  only  one  oxidant  channel 
and  adjacent  rib  are  included.  A  computational  model  built  in  this 
way  can  not  only  deal  with  fuel  cells  with  oxidant  channels  hav¬ 
ing  straight,  serpentine,  and  Z  patterns  but  can  also  be  applied  to 
large-scale  fuel  cells  while  using  fewer  resources  and  less  time.  In 
the  following  section,  the  basic  assumptions  used  in  this  study  are 
elucidated. 

2.2.  Basic  assumptions 

The  dimensions  considered  by  this  model  can  be  illustrated  with 
Fig.  1(a)  and  (b).  In  Fig.  1(a),  only  the  dimension  along  the  direction 
of  oxidant  flow  is  adopted;  thus,  this  model  can  be  regarded  as  one¬ 
dimensional  in  the  x-y  plane.  In  Fig.  1(b),  only  the  dimension  along 
the  direction  normal  to  active  areas  is  adopted;  hence,  this  model 
can  be  regarded  as  one-dimensional  in  the  y-z  plane.  This  model 
consisting  of  one  dimension  (along  the  direction  of  oxidant  flows) 
plus  one  dimension  (along  the  z-direction)  forms  the  main  frame¬ 
work  for  developing  the  following  mathematical  formulations.  The 
basic  assumptions  were  as  follows: 

1.  Overpotentials  within  the  anode  catalyst  layer  (ACL)  are 

ignored,  and  only  those  within  the  CCL  are  considered. 
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Fig.  1.  (a)  Fuel-cell  computational  models  in  the  x-y  plane,  (b)  Fuel-cell  computa¬ 
tional  models  in  the  y-z  plane. 


2.  According  to  assumption  1,  the  electrochemistry  and  trans¬ 
port  phenomena  at  the  anode  can  be  ignored;  only  those  at 
the  cathode  are  considered. 

3.  Flow  fields  at  the  cathode  are  steady  and  laminar. 

4.  Gas-phase  mixtures  at  the  cathode  are  regarded  as  ideal  gases. 

5.  Non-uniform  distributions  of  oxidant  flows  are  ignored;  thus, 
the  oxidant  can  be  regarded  as  feeding  each  channel  uniformly. 

6.  For  mass  transfer  of  each  gas-phase  species,  only  convection 
along  the  streamwise  directions  within  the  channels,  and  only 


diffusion  along  the  z-direction  within  other  domains  are  con¬ 
sidered. 

7.  Only  gas-phase  pressure  variations  along  streamwise  direc¬ 
tions  within  the  channels  are  considered,  and  gas-phase 
pressure  gradients  along  the  z-direction  within  other  domains 
are  regarded  as  zero. 

8.  The  temperatures  on  the  rib  boundaries  are  specified,  and 
within  other  domains,  only  heat  conduction  along  the  z- 
direction  is  considered. 

9.  Only  current  transfer  along  the  z-direction  is  considered,  and 
current  densities  are  assumed  to  be  uniform  in  the  z  direction. 

10.  According  to  assumption  9,  both  electronic-phase  and  ionic- 
phase  potential  variations  only  occur  in  the  z-direction. 

11.  Only  liquid-water  transport  along  the  streamwise  direc¬ 
tions  within  channels  and  along  the  z-direction  within  other 
domains  is  considered. 

1 2.  Property  variations  along  the  direction  normal  to  oxidant  flows 
in  the  x-y  plane  are  ignored;  thus,  properties  are  assumed  to 
be  constant  in  this  direction. 

13.  Property  variations  in  the  z-direction  on  PEM  boundaries  in  the 
vicinity  of  the  anode  are  ignored;  thus,  on  these  boundaries,  all 
property  gradients  in  the  z  direction  are  zero. 

2.3.  Mathematical  formulation 


2.3  A.  Mass  conservation  of  each  gas-phase  species 
di  dz\gdz 


+  S‘ 


(1) 


Eqs.  (1)  and  (2)  describe  mass  conservation  of  each  gas-phase 
species  within  oxidant  channels  and  other  domains,  respectively.  In 
these  equations,  l  and  z  represent  coordinates  along  the  streamwise 
direction  (the  direction  of  oxidant  flows)  and  z-direction,  respec¬ 
tively;  i  is  the  gas-phase  species  of  oxygen,  nitrogen,  or  water  vapor; 
Ufj  is  the  gas-phase  molar  flux  along  the  streamwise  direction;  D'g 
is  the  gas-phase  diffusivity;  c,-  is  the  gas-phase  molar  density;  and 
Slc  is  the  source  term  of  each  gas-phase  species. 

Eqs.  ( 1 )  and  (2)  can  be  simplified  by  expressing  cz  with  gas-phase 
partial  pressure  as  follows: 


where  Pg  is  the  gas-phase  partial  pressure,  R  is  the  ideal-gas  con¬ 
stant,  and  T  is  the  absolute  temperature.  Substituting  Eq.  (3)  into 
Eqs.  (1)  and  (2)  gives 
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In  this  study,  Eqs.  (4)  and  (5)  are  adopted  to  describe  the  mass  con¬ 
servation  of  each  gas-phase  species  within  oxidant  channels  and 
other  domains,  respectively.  Using  the  Bruggeman  correction,  Dlg 
can  be  expressed  as  [29]  follows: 


)f  _  £T(1  -s)r(l  -x}g) 
jj  ±  i 


(6) 


where  s  is  the  porosity,  s  is  the  liquid  saturation,  r  is  the  tortuosity, 
xlg  is  the  gas-phase  molar  fraction,  and  D|  is  the  binary  diffusivity  of 
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both  the  i-th  and  j-th  gas-phase  species,  which  can  be  determined  decided  by  the  ratio  of  -{dcH2oldz)c  to  (dcH2oldz)a  -  ( dcH2o/dz)c . 
as  follows  [30]:  The  terms  ( dcu2oldz)a  and  (dcu2oldz)c  can  be  expressed  as  follows: 


ug  ~ 
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where  a  and  b  are  empirical  constants  and  P2J-  is  the  sum  of  par¬ 
tial  pressures  of  both  the  i-th  and  j-th  gas-phase  species.  Pd  and 
Pcj  are  the  critical  pressures  for  the  i-th  and  j-th  gas-phase  species, 
respectively;  Td  and  TCJ  are  the  critical  temperatures  for  the  i-th  and 
j-th  gas-phase  species,  respectively;  and  M2  and  Mj  are  the  molec¬ 
ular  weights  for  the  i-th  and  j-th  gas-phase  species,  respectively.  In 
addition,  Slc  can  be  expressed  as  follows: 
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Eq.  (8)  defines  the  source  term  of  oxygen,  S°2,  where  jj  is  the  vol¬ 
umetric  transfer  current  density,  and  F  is  the  Faraday  constant.  It 
can  be  determined  from  Eq.  (8)  that  within  the  CCL,  S°2  represents 
the  oxygen  consumption  rate  caused  by  electrochemical  reactions, 
while  within  other  domains,  this  term  is  zero.  For  5^2 ,  because  the 
nitrogen  does  not  participate  in  electrochemical  reactions,  it  is  nei¬ 
ther  consumed  nor  generated,  which  leads  to  S^2  being  zero  within 
all  domains.  Eq.  (10)  describes  the  source  term  of  water  vapor,  5^2°f 
where  Ss  is  the  source  term  of  liquid  water  and  thus  -Ss  is  the  sink 
term  of  water  vapor,  and  MH2o  is  the  molecular  weight  of  water. 
yj^Q  is  the  correction  factor,  which  indicates  the  ratio  of  transport 
into  the  cathode  to  total  transport  for  water  vapor  produced  in  elec¬ 
trochemical  reactions.  In  addition,  X  is  the  water  content,  and  iz  is 
the  current  density  in  the  z  direction. 

It  can  be  determined  from  Eq.  (10)  that  this  term  -Ss/M h2o 
exists  within  all  domains,  depending  on  the  phase  changes  between 
water  vapor  and  liquid  water.  Within  the  CCL,  water  vapor  is  pro¬ 
duced  at  the  rate  of  jr/2F  according  to  electrochemical  reactions 
in  which  the  transport  of  water  vapor  toward  the  anode  and  cath¬ 
ode  is  related  to  the  gradient  of  water-vapor  molar  density  at  the 
anode  and  cathode,  respectively.  In  order  to  evaluate  the  ratio  of 
water-vapor  transport  toward  the  cathode  gas  channel,  y^Q  can 
be  defined  as  follows: 


_1  d  1 

R  dz  1 

{  —  )\ 

o 


for 

for 


\1 

d  i 

/PH2°X1 

R 

dz 

l  T  ) 

1 

d  1 

(  Pg2°  \  " 

R 

dz  1 

K~)\ 

<0 

>  0 


(12) 


'  V%&-r5&)/Krc* 

dpEM  +  ( dAa  +  dca.)/  2 

0 


for 

for 


(PgH&-p5qyRT^ 

dpEM  +  ( dAa  +  dca)l  2 

dpEM  +  (dACL  +  dCa)l  2 


(13) 


In  Eq.  (12),  {dcu2oldz)c  can  be  rewritten  as  [(1  /R)d(P^2°/T)/dzj  c 
according  to  Cu2o  =  Pg2°  l{RT).  When  [(l/F)d(P|I20/T)/dz]c  <  0, 
i.e.,  water  vapor  diffuses  from  the  CCL  toward  the  cathode, 
{dcH20ldz)c  is  equal  to  [(l/fl)d(P|h°/r)/dz]c.  However,  when 
[(l/P)d(P|I2°/T)/dz]c  >  0,  i.e.,  water  vapor  diffuses  from  the  cathode 
toward  the  CCL,  {dcu2oldz)c  is  forced  to  be  zero,  which  assumes  that 
it  is  impossible  for  water  vapor  to  diffuse  from  the  cathode  to  the 
CCL.  In  Eq.  (13),  P^2C°L  and  P^o.  are  water- vapor  partial  pressures 
within  the  CCL  and  ACL,  respectively,  where  Pg^L  can  be  obtained 
from  iterative  calculations,  while  j°L  needs  to  be  evaluated  in  an 
approximate  manner  because  the  electrochemistry  and  transport 
phenomena  at  the  anode  are  ignored.  In  addition,  TCcl  is  the  tem¬ 
perature  within  the  CCL  approximated  as  that  within  the  ACL,  and 
dpEM >  dAcL,  and  dCa  are  the  thicknesses  of  the  PEM,  ACL  and  CCL, 
respectively.  It  can  be  determined  from  Eq.  (13)  that  (dcH2o/dz)a 
can  be  determined  by  the  relation  of  {Pg2c^L  -  P^°l)/PTccl  divided 
by  dpEM  +  (dACL  +  dca )/2,  where  (Pg}°L  -  PfA°L  )/RTccl  represents  the 
differences  between  the  water-vapor  molar  density  within  the  CCL 
and  the  ACL.  When  this  relation  is  larger  than  zero,  i.e.,  water 
vapor  diffuses  from  the  CCL  toward  the  anode,  (dcH2o/dz)a  can  be 
obtained  by  this  relation.  However,  when  this  relation  is  less  than 
zero,  {dcu2oldz)a  is  forced  to  be  zero,  which  assumes  that  water- 
vapor  diffusion  from  the  anode  to  the  CCL  is  impossible.  In  order  to 
evaluate  P^ j°L,  an  approximate  method  is  adopted,  as  follows. 

First,  the  water-vapor  partial  pressure  at  the  anode  inlet,  P^1  ? °a, 
can  be  calculated  as  follows: 


PXa=RH?nPsaia 


(14) 


where  RHlan  and  P™f  are  relative  humidity  and  saturation  pressure 
at  the  anode  inlet,  respectively.  According  to  ideal-gas  laws,  the 
hydrogen  molar  flow  rate  at  the  anode  inlet,  N|^a,  can  be  expressed 
as  follows: 


rptotal  _  pH^O  )Qin 
j^in,a  _  K  g,in,a  g,in,aJ^a 

H2  “  RTlan 


(15) 


vEC  - 
Kh,o  - 


0  for(dcH2o/dz)a  =  0  and  (dcH2o/dz)c  =  0 


-{dcu2oldz)c 
{dcH2oldz)a  -  {dcu2oldz)c 


otherwise 


(11) 


In  Eq.  (11),  ( dcu2oldz)a  and  ( dcu2oldz)c  represent  the  gradients  in 
the  z-direction  of  water- vapor  molar  density  at  the  anode  and  cath¬ 
ode,  respectively.  When  (dcu2oldz)a  and  ( dcH2o/dz)c  are  both  zero, 
this  represents  the  fact  that  water  vapor  produced  by  electrochem¬ 
ical  reactions  neither  diffuses  from  the  CCL  toward  the  anode  gas 
channel  nor  toward  the  cathode  channel,  thus  leading  to  yEE0  being 

zero.  However,  when  one  of  these  two  terms  is  not  zero,  it  means 
that  water  vapor  produced  by  electrochemical  reactions  may  dif¬ 
fuse  from  the  CCL  toward  the  anode  or  toward  the  cathode,  where 
diffusion  fluxes  toward  the  anode  and  cathode  are  proportional  to 
( dc\\2oldz)a  and  ~(dcu2oldz)c1  respectively.  Therefore,  yEE0  can  be 


In  Eq.  (15),  the  gas-phase  total  pressure  at  the  anode  inlet,  Pgf°la 

is  approximated  as  that  at  the  anode  outlet,  P|°^t  a,  and  Tan  is  the 
temperature  at  the  anode  inlet.  Therefore,  substituting  the  known 
values  N[^a,  Pg°f“la,  T ™  and  P^f°a  into  Eq.  (15)  gives  the  gas-phase 

volumetric  flow  rate  at  the  anode  inlet,  Q^n.  After  obtaining  Q "h 
according  to  ideal-gas  laws,  water-vapor  molar  flow  rates  at  the 
anode  inlet  and  outlet,  i.e.,  N^q  and  can  be  expressed  as 

follows: 


PH ?°  Qj1 

ATout,a  ~  \jin,a  _  g,in,a^a 
iVH20  ~  JVH20  -  RJin 


(16) 


In  Eq.  (16),  N^a  is  approximated  as  N|^q,  which  assumes  that 
water-vapor  molar  flow  rates  within  the  fuel  channel  are  constant. 
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The  hydrogen  molar  flow  rate  at  the  anode  outlet,  N°"t,a,  can  be 
expressed  as  follows: 


i 

2 F 


(17) 


where  I  is  the  output  current,  and  thus  //2F  stands  for  the  total 
consumption  rate  of  hydrogen  in  power  generation. 

With  Eqs.  (15)-(17),  the  water-vapor  partial  pressure  at  the 
anode  outlet,  can  calculated  as  follows: 


p  total  iyOut,a 
H20  =  1  g,out,a^H20 

1  g,out,a  ArOut,a  ,  *  rOut,a 
iVH2  +JVH20 


(18) 


After  obtaining  ? °a  and  Pg20°t  a,  Pg  can  be  evaluated  using 
linear  interpolation  between  Pg?°a  and  as  follows: 


pH20  _ 
1  g,ACL  ~ 


P"]n°a  +  ^(P"fut,a-PXa)f0r  Parallel  fl°W 

Pg,iut,a+^(Pgjn,a-pg,2out,Jfm  counter  flow 


(19) 


pg  can  be  expressed  as  follows: 


N 

5>,M, 


Pg  ~ 


VgJ 


(24) 


where  /Mj  is  the  sum  of  mass  fluxes  of  each  gas-phase  species. 

i 

fg  can  be  expressed  as  [31  ]  follows: 


fg  ~ 


55+4,-5ex[,(^f) 


Rpgkrg 


(25) 


where  Reg  is  the  gas-phase  Reynolds  number  and  krg  is  the  gas- 
phase  relative  permeability.  Reg  can  be  defined  as  follows: 


Peg  - 


Pg  Vg./fl 

P'g 


(26) 


where  iig  is  the  gas-phase  dynamic  viscosity  and  can  be  determined 
as  follows  [30]: 


where  Lp  is  the  partial  length  from  the  inlet  to  the  corresponding 
computational  grids  of  the  oxidant  channel  and  Lt  is  the  total  length 
of  the  oxidant  channel.  When  the  flows  of  fuel  and  oxidant  are 
parallel,  the  first  term  in  Eq.  (19)  is  used  to  evaluate  P^ 2°L,  and 
when  flows  of  fuel  and  oxidant  are  counter  to  one  another,  the 
second  term  in  Eq.  (19)  is  used. 

In  addition,  it  can  be  determined  from  Eq.  (10)  that  within 
the  CCL  and  PEM,  water-vapor  transport  driven  by  proton  drags 
d[(- 2.5A,/22)(zz/F)]/dz  also  needs  to  be  considered  [1].  This  term 
describes  how  when  protons  move  with  the  rate  of  zz/F  from  the 
ACL  to  the  CCL,  water  vapor  dragged  by  these  ions  can  also  move 
together  with  them,  where  the  number  of  water  molecules  carried 
by  one  proton  is  equal  to  2.5A/22.  X  is  further  defined  below. 


2.3.2.  Gas-phase  momentum  conservation 


dPg 

dl 


’<L,g  PtV, g.l 
SI  2 


+  pgg  cos  e,/g 


(20) 


Eq.  (20)  describes  gas-phase  momentum  conservation  within  the 
oxidant  channel,  where  Pg  is  the  gas-phase  pressure,  l  is  the  coor¬ 
dinate  along  the  streamwise  direction,  fg  is  the  gas-phase  friction 
factor,  and  D  is  the  channel  hydraulic  diameter.  pg  is  the  gas-phase 
density,  Vgi  is  the  gas-phase  velocity  along  the  streamwise  direc¬ 
tion,  I<Ltg  is  the  gas-phase  local  resistance,  and  81  is  the  local  channel 
length.  In  addition,  g  is  the  gravitational  acceleration,  and  0i/g  is  the 
included  angle  of  oxidant  flows  relative  to  gravity,  where  Oi/g  -  0° 
and  Qi/g  =  180°  represent  the  fact  that  oxidant  flows  are  parallel 
and  counter  to  gravity,  respectively,  while  0ijg  =  90°  represents 
the  fact  that  oxidant  flows  cannot  be  affected  by  gravity.  D  can  be 
expressed  as  follows: 


2wcdc 
wc  +  dc 


(21) 


where  wc  is  the  channel  width  and  dc  is  the  channel  depth. 
Vgi  can  be  determined  as  follows: 


XlgPilg 


.y  1  ( i  Mi 


-1/2 


1  + 


(27) 


(28) 


where  p}g  and  pdg  are  the  dynamic  viscosities  of  the  z-th  and  j-th 
gas-phase  species,  respectively,  which  can  be  calculated  according 
to  Sutherland’s  law  as  follows: 


pilg  = 


AjT1-5 

VTT 


(29) 


where  Az  and  B/  both  are  empirical  constants,  and  T  is  the  absolute 
temperature. 

Krg  can  be  expressed  as  [32,33]  follows: 


Krg  —  (1  —  s)3 


(30) 


It  can  be  determined  from  Eq.  (30)  that  when  s  is  zero  (i.e.,  there 
is  an  absence  of  liquid  water  within  the  channel),  krg  is  one,  while 
when  s  is  not  zero  (i.e.,  the  presence  of  liquid  water  within  the 
channel),  Krg  decreases  with  increasing  s.  Therefore,  combining  Eqs. 
(25)  and  (30)  shows  that  fg  can  increase  with  increasing  s,  which 
means  that  stronger  flooding  can  lead  to  greater  gas-phase  flow 
resistance  produced  within  the  channel. 

I<Ltg  can  be  expressed  as  [34]  follows: 


f  0  for  straight  channel  sections 
^L  g  \  Jjfs  f°r  channel  bends 

It  can  be  determined  from  Eq.  (31)  that  \<itg  is  zero  for  straight 
channel  sections,  while  I<Ltg  can  be  written  as  the  product  of  Le/D  and 
fg  for  channel  bends,  where  Le/D  is  the  dimensionless  equivalent 
length  depending  on  the  relative  radius  of  channel  bends. 


vg.i  = 


tigjRT 

~pT 


(22) 


where  ng  [  is  the  gas-phase  molar  flux  along  the  steamwise  direction 
and  can  be  expressed  with  the  sum  of  molar  fluxes  of  each  gas- 
phase  species  as  follows: 


N 

ns,i = J2nt’1 


(23) 


2.3.3.  Energy  conservation 

E(%S)+ST.0  (32) 

Eq.  (32)  describes  heat  conduction  in  the  z-direction  within  the 
fuel  cell,  where  keff  is  the  effective  thermal  conductivity  and  ST  is 
the  heat  source.  keffC an  be  expressed  as  [29]  follows: 

K eff  —  —  2  ks  + 


s/ (2  ks  +  kf 

mix  )  +  l-e/(3  ks) 


(33) 
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where  ks  and  kfmix  are  the  thermal  conductivities  of  solids  and 
fluids,  respectively,  and  kfmix  can  be  defined  as  follows: 

kfmix  =  s^l  +  (1  —  s)kg  (34) 

where  /q  and  kg  are  the  thermal  conductivities  of  liquid  phases  and 
gas  phases,  respectively,  and  kg  can  be  defined  as  [30]  follows: 


within  the  CCL,  respectively.  Pg2  can  be  obtained  from  iterative 
computations,  but  P|b  needs  to  be  evaluated  using  Eq.  (19)  and 
PgACL  =  Pg°out,a  ~  pgACL •  open-circuit  voltage  within  the  ACL, 
VqC,  is  assumed  to  be  zero. 

Overpotentials  can  be  determined  according  to  the  following 
equations  [6,31]: 


ka  — 


N 

=£- 


xgkg 

E/g0! 


Mj 


-1/2 


1  + 


(35) 


(36) 


where  klg  and  kfg  are  the  thermal  conductivities  of  the  i-th  and  j-th 
gas-phase  species,  respectively,  which  can  be  determined  accord¬ 
ing  to  the  Eucken  formula  [30]  as  follows: 

l<lg  =  ^ Cjj  g  +  -  V's 

where  Clp  g  is  the  constant-pressure  specific  heat  of  the  i-th  gas- 
phase  species  and  can  be  determined  according  to  the  Jannaf 
method  as  follows: 

CP,g  =  (“'l  +  a27  +  Q3T2  +  ^  +  Q5T4)  08) 


where  a\  -  a*5  are  empirical  constants  and  T  is  the  absolute  tem¬ 
perature. 

ST  can  be  expressed  as  follows: 


ST 


-iz~!r-  within  the  rib  and  CGDL 
dz 

<  ^hin  the  CCL 


within  the  PEM 


(39) 


where  4>s  and  are  electronic-phase  and  ionic-phase  potentials, 
respectively,  and  p  is  the  overpotential.  It  can  be  determined  from 
Eq.  (39)  that  heat  sources  within  both  the  rib  and  the  CGDL  are 
Joule  heating  equal  to  the  product  of  iz  and  -d</>s/dz,  which  repre¬ 
sents  the  conversion  from  electric  energy  to  heat  energy  caused  by 
electric-phase  resistivity.  The  heat  source  within  the  CCL  includes 
not  only  the  Joule  heating  but  also  the  irreversible  heat  source, 
jTr /,  which  represents  the  waste  heat  produced  in  the  process  of 
conversion  from  chemical  energy  to  electric  energy.  Heat  sources 
within  the  PEM  are  the  Joule  heating  equal  to  the  product  of  iz  and 
-d0j/dz,  which  represents  the  conversion  from  electric  energy  to 
heat  energy  caused  by  ionic  resistivity. 

Interfacial  temperature  changes  resulting  from  contact  thermal 
resistance  between  the  rib  and  CGDL  can  be  calculated  as  follows: 


ST  =  R’)blCDL{izS(prsiblCDL) 


(40) 


where  ST  is  the  interfacial  temperature  change,  R'fGDL  is  the  con¬ 
tact  thermal  resistance  of  the  rib/CGDL  interfaces,  and  <S0+GDL 
is  the  interfacial  change  of  the  electric-phase  electric  potentials 
caused  by  RjblGDL. 


2.3.4.  Open-circuit  voltages  and  polarization  losses 

Open-circuit  voltages  can  be  expressed  as  [31  ]  follows: 

VGC  =  1.482  -  0.000845T  +  0.000043 IT ln[P^2(p02)°-5]  (41) 

Eq.  (41)  defines  the  open-circuit  voltage  within  the  CCL,  V^c, 
where  T  is  the  absolute  temperature  and  P|h  and  P° 2  are  the  hydro¬ 
gen  partial  pressure  within  the  ACL  and  the  oxygen  partial  pressure 


iz 


(1  -s)ir0ef 


X 


-  exp 


( 


-aFr]c  \ 
RT  ) 


(42) 


where  iz  is  the  current  density  in  the  z-direction,  (1  -  s)  represents 
the  ratio  of  effective  catalyst  surface  area  to  total  catalyst  surface 
area,  is  the  reference  transfer  current  density,  Pg2  is  the  oxy¬ 
gen  partial  pressure,  Pr^  is  the  reference  pressure,  and  Ec  is  the 
activation  energy  for  oxygen  reduction  on  Pt.  Tref  is  the  reference 
temperature,  a  is  the  transfer  coefficient,  and  r\c  is  the  cathode 
overpotential.  iz  can  be  calculated  with  Eq.  (42),  and  iz/dca  gives  jj 
because  jT  can  be  approximated  as  the  gradient  of  iz  in  the  z  direc¬ 
tion,  i.e.,jY  =  diz\dz.  During  iterative  computations,  in  order  to  solve 
rjc ,  another  expression  equivalent  to  Eq.  (42)  needs  to  be  used: 


RT  . 


i+(i+1 


(43) 


T]c  determined  with  Eq.  (43)  is  positive,  where  i0)Z  is  the  surface 
transfer  current  density  and  can  be  expressed  as  follows: 

/  \  0.5 


*  £ l 


<0,Z  =  (1-S)CM^ 


(44) 


Ionic-phase  ohmic  losses,  8c/)f,  can  be  calculated  as  follows: 

8<t>f  =  l^~  (45) 

Of 

8(/)f  determined  with  Eq.  (45)  is  positive,  where  df  is  the  effective 
thickness  equal  to  dPEM  +  ( dAa  +  dccr)/2. 6/  is  the  ionic  conductivity 
and  is  defined  as  the  average  of  a^CL  and  o^a,  i.e.,  oy  =  (o^a  + 
<jjCL)/ 2.  Of  at  the  anode  or  cathode  can  be  expressed  as  [1  ]  follows: 


Of  =  (0.5139A 


0.326)  exp 


1268 


k  = 


0.043  +  17.18a  -  39.85a2  +  36.0a3  0  <  a  <  1 
14  +  1.4(a-l)  1  <  a  <  3 

16.8  a  >  3 


p_t °_ 

Psat 


(46) 

(47) 

(48) 


p  _  „-5800.2206T-1+1.3914993-0.048640239T+0.41764768x10-4T2-0.14452093x10-7T3+6.5459673  In 7 

1  sat  —  C 

(49) 

It  can  be  determined  from  Eq.  (46)  that  Of  is  a  function  of  water  con¬ 
tent,  A  and  absolute  temperature,  T.  Eq.  (47)  shows  that  A  depends 
on  water  activity,  a,  which  is  related  to  water-vapor  partial  pres¬ 
sure,  P|b0  and  saturation  pressure,  Psat,  as  shown  in  Eq.  (48),  where 
Psat  can  be  calculated  using  Eq.  (49)  [31].  When  Eqs.  (46)-(49) 
are  used  to  determine  ojCL,  all  of  the  necessary  properties  can  be 
obtained  from  iterative  computations.  However,  for  ionic  conduc¬ 
tivity  within  the  ACL,  Pgj^L  from  Eq.  (19)  and  T,  approximated  as 
the  temperature  within  the  CCL,  are  substituted  into  Eqs.  (46)-(49) 
to  evaluate  OjCL. 
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Electronic-phase  ohmic  losses,  8(ps  can  be  expressed  as  follows: 


8<ps  =  izRs,z 


(50) 


8(j)s  determined  from  Eq.  (50)  is  positive,  where  RSyZ  is  the  effective 
electric  resistivity  in  the  z-direction  of  the  bipolar  plate  or  GDL,  or 
the  effective  contact  electric  resistance  in  the  z-direction  between 
the  rib  and  GDL,  which  can  be  expressed  as  [31  ]  follows: 

BP  fwL+wc  dBp\  (Wi+WC)WC 

=  l^^dc  +  T )  Rbp’z  + - 4dgp - Rsp’xy  (51 } 


oCDL  w  D  ,  (Wl+WC)WCd 

Ks,z  =  aGDLKGDL,z  H - - KGDL,xy 

ougdl 

dCR  _  WL  +  Wc  DriblGDL 
s’z  “  WL  CR 


(52) 

(53) 


Eqs.  (51 )  and  (52)  define  the  effective  electric  resistivity  of  the  bipo¬ 
lar  plate  and  GDL,  Rpp  and  R^L,  respectively,  and  Eq.  (53)  defines 
the  effective  contact  electric  resistance  of  rib/GDL  interfaces,  R£R.  In 
these  equations,  wL  is  the  rib  width  between  gas  channels  or  coolant 
channels,  wc  is  the  width  of  gas  channels  or  coolant  channels,  dG  is 
the  depth  of  gas  channels  or  coolant  channels,  dBp  is  the  thickness  of 
the  bipolar  plate  minus  the  sum  of  the  depth  of  gas  channels  plus 
that  of  coolant  channels,  and  dGDL  is  the  GDL  thickness.  RBp,z  and 
RBPyXy  are  the  electric  resistivity  of  the  bipolar  plate  in  the  z  and  x-y 
directions,  and  r!^gdl  is  the  contact  electric  resistance  of  rib/GDL 
interfaces. 

For  the  bipolar  plate  with  one  surface  containing  gas  channels 
and  another  one  containing  coolant  channels,  Rpp  needs  to  be  eval¬ 
uated  for  the  gas  and  coolant  side,  respectively,  while  for  the  bipolar 
plate  with  only  one  surface  containing  gas  channels,  Rpp  needs  to 
be  evaluated  only  for  the  gas  side.  To  calculate  the  electronic-phase 
ohmic  losses  of  a  single  cell,  Eqs.  (51)-(53)  are  first  used  to  evalu¬ 
ate  Rrp,  R^l  and  R^r  at  the  anode  and  cathode,  respectively.  Then, 
Eq.  (50)  is  applied  to  both  the  anode  and  cathode  to  obtain  8(ps- 
However,  in  this  study,  because  only  current  densities  at  the  cath¬ 
ode  side  are  available,  electronic-phase  ohmic  losses  at  the  anode 
or  coolant  side  are  evaluated  by  substituting  the  average  current 
densities  into  Eq.  (50).  In  addition,  it  can  be  determined  from  Eqs. 
(51 )— (53)  that  effective  electric  resistivity  depends  not  only  on  the 
properties  of  the  component  itself,  but  also  the  effects  of  geome¬ 
try  parameters  on  current  transfer.  For  example,  the  bipolar  plate 
with  wider  and  deeper  channels  can  yield  greater  Rpp,  and  when 
the  GDL  is  thicker  or  the  gas  channel  is  wider,  R^ZL  becomes  larger. 
RribiGDL  can  jncrease  wjti1  increasing  ratios  of  the  channel  width  to 
rib  width,  wc/wL. 


2.3.5.  Computations  of  local  cell  voltages 

In  order  to  compute  local  cell  voltages,  this  study  presents 
an  equation  capable  of  describing  how  local  cell  voltages  are 
related  to  the  open-circuit  voltage  and  various  polarization  losses. 
Fig.  2  illustrates  the  variations  of  electronic-phase  and  ionic-phase 
potentials  within  the  cell.  It  can  be  determined  from  Fig.  2  that 
the  electronic-phase  potential  at  the  anode  boundary  is  and 
the  electronic-phase  ohmic  loss  at  the  anode  is  8<pf,  thus  lead¬ 
ing  to  the  electronic-phase  potential  within  the  ACL,  0^a,  which 
is  equal  to  <p£  -  8q bf.  The  correlation  between  the  open-circuit 
voltage,  overpotential  and  both  electronic-phase  and  ionic-phase 
potentials  within  the  ACL  can  be  expressed  as  follows: 

-riA  =  $CL  -  <PfCL  ~  Voc  (54) 


PEM 


u 


Fig.  2.  Variations  of  electronic-phase  and  ionic-phase  potentials  within  the  fuel  cell. 


of  ionic-phase  potentials  between  the  ACL  and  CCL  can  be  expressed 
as  follows: 

tfCL  =  tfCL- 8. 4>f  (55) 

where  4>jCL  is  the  ionic-phase  potential  within  the  CCL  and  8</>f  is 
the  ionic-phase  ohmic  loss  between  the  ACL  and  CCL.  The  corre¬ 
lation  between  the  open-circuit  voltage,  overpotential  and  both 
electronic-phase  and  ionic-phase  potentials  within  the  CCL  can  also 
be  expressed  in  a  manner  similar  to  Eq.  (54): 

~nc  =  <t>sCL  -  4>jCL  -  Vo,  (56) 

where  r)C  is  the  cathode  overpotential  and  (p^CL  is  the  ionic-phase 
potential  within  the  CCL.  Combining  Eqs.  (54)-(56)  yields  (pGCL‘. 

<PsCL  =  V0Cc  +  #  -  S#  -  %  -  »/c  (57) 

After  obtaining  0fa,  the  electronic-phase  potential  at  the  cath¬ 
ode  boundary  can  be  written  as  follows: 

0sC  =  (PsCL  -  WCS  (58) 

where  is  the  electronic-phase  potential  at  the  cathode  boundary 
and  8(ps  is  the  electronic-phase  ohmic  loss  at  the  cathode.  There¬ 
fore,  the  local  cell  voltage,  </>f  -  (pf,  can  be  expressed  as  follows: 

<pcs  -  =  v£c  -  rtc  -  8<ps  -  8<pf  (59) 

where  8(ps  is  the  sum  of  the  electronic-phase  ohmic  losses  at  the 
anode  and  cathode,  i.e.,  8(ps  =  8<p£  +  8<p%.  It  can  be  determined  from 
Eq.  (59)  that  the  local  cell  voltage  is  equal  to  the  open-circuit  volt¬ 
age  minus  all  polarization  losses.  According  to  Eq.  (59),  the  local  cell 
voltage  can  be  obtained  by  evaluating  Vfc,  pc ,  &4>s  and  8<pf ,  respec¬ 
tively.  Compared  with  solving  potential  equations  with  the  Laplace 
form,  this  study  uses  simple  algebraic  operations  in  place  of  com¬ 
plicated  differential  equations,  thus  numerical  computations  can 
be  further  simplified. 


where  pA  is  the  anode  overpotential,  assumed  to  be  zero,  0^CL  is  the 
ionic-phase  potential  within  the  ACL,  and  v£c  is  the  open-circuit 
voltage  within  the  ACL,  and  its  value  is  zero.  Therefore,  it  can  be 
determined  from  Eq.  (54)  that  0^CL  is  equal  to  (pjCL.  The  correlation 


2.3.6.  Liquid-water  transport  and  production 

In  order  to  simulate  two-phase  flows  within  the  fuel  cell,  the 
two-phase  Darcy’s  law  [32],  originally  only  applied  to  porous 
media,  is  further  developed  into  mathematical  formulations  that 
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can  describe  liquid-water  transport  and  production  not  only  within 
porous  media  but  also  within  the  oxidant  channel.  By  solving 
related  governing  equations,  liquid  saturation  is  obtained  so  that 
the  effects  of  flooding  on  the  local  phenomena  and  the  overall  per¬ 
formance  of  fuel  cells  can  be  evaluated.  According  to  Darcy’s  law 
for  two-phase  flow,  liquid-phase  momentum  conservation  within 
the  oxidant  channel  can  be  expressed  as  follows: 


dp, 

dl 


viPiVu 

Kpkfi 


+  pig  cos  9t /g 


(60) 


where  P/  is  the  liquid-phase  pressure,  /  is  the  coordinate  along  the 
streamwise  direction,  vt  is  the  liquid-phase  kinematic  viscosity,  and 
Kfj  is  the  equivalent  absolute  permeability  of  the  channel.  p/V/  /  is 
the  liquid-phase  mass  flux  along  the  streamwise  direction,  where 
Pi  is  the  liquid-phase  density,  and  Vitl  is  the  liquid-phase  velocity 
along  the  streamwise  direction.  kri  is  the  liquid-phase  relative  per¬ 
meability,  which  can  be  determined  according  to  kri  =  s3  [33  ],  and  g/ 
is  the  component  in  the  l  direction  of  gravitational  acceleration.  P/ 
can  be  expressed  with  the  gas-phase  pressure,  Pg,  and  the  capillary 
pressure,  Pc,  as  [32]: 


Pi  =  Pg-Pc 


where  Pc  can  be  determined  as  follows  [35]: 


Pc 


a  cos  0C 


1/2 

Ks) 


(61) 


(62) 


For  hydrophobic  porous  media  of  0C  <  90° : 


ViPM'Z 

I<p 


-s3  p^  cos  ez/g 


(67) 


where  I<p  is  the  absolute  permeability  of  porous  media,  p/V/z  is 
the  liquid-phase  mass  flux  in  the  z-direction,  where  Viz  is  the 
liquid-phase  velocity  in  the  z-direction,  and  6z/g  is  the  included 
angle  of  liquid  flows  relative  to  the  gravity  within  porous  media. 
( Oz/g  =  0°  and  0zig  =  1 80°  represent  the  fact  that  liquid  flows  are  par¬ 
allel  and  counter  to  gravity,  respectively,  while  0z/g  =  90°  represents 
the  fact  that  liquid  flows  cannot  be  affected  by  gravity.) 

In  order  to  calculate  Kjj  in  Eqs.  (64)  and  (65),  first,  Eq.  (60)  can 
be  rewritten  in  a  manner  similar  to  Eq.  (20)  as  follows: 


dPi 

dl 


f,  pK 

D  2 


Kl,1  pK 

SI  2 


+  pig  cos  e[/g 


(68) 


where  the  friction  factor,//  can  be  expressed  as  follows: 


fi= 


55  +41.5  exp 


\wc/dcJ 


Reikri 


(69) 


The  liquid-phase  Reynolds  number,  Reh  can  be  expressed  as  fol¬ 
lows: 


Re  i  = 


PiVijD 

\ii 


(70) 


The  liquid-phase  local  resistance,  I<Lh  can  be  expressed  as  follows: 


i(s)  = 


f  1.417(1  -s)- 2.120(1  -s)2  + 1.263(1  -s)3, 
\  1.417s  -  2.120s2  +  1.263s3, 


if  0c  <  90° 
if  0c  >  90° 


(63) 


In  Eq.  (62),  o  is  the  surface  tension,  0C  is  the  contact  angle,  s 
is  the  porosity,  and  /(s)  is  the  Leverett  J-function,  as  shown  in  Eq. 
(63).  Combining  Eqs.  (60)-(63)  and  rearranging  yields  governing 
equations  with  the  variable  s. 

For  the  hydrophilic  channel  of  0C  <  90° : 


For  the  hydrophobic  channel  of  0c  >  90° : 


{0  for  channel  straight  sections 
for  channel  bends 

Substituting  Eqs.  (69)-(71)  into  Eq.  (68)  gives 


(71) 


dPi  _  [55  +  41 .5  exp(-3.4/(wc/dc))](l  +  rI<LD/8l)  V/P/VM 

dl  2D2  1^ 


+  P/g  cos  %  (72) 

Then,  by  comparing  equivalent  Eqs.  (60)  and  (72),  K*  can  be  deter¬ 
mined  as  follows: 


lp  [55  +  41.5  exp(-3.4/( wc/dc))]  (1  +  (f*tD/5/)) 

Jo  for  channel  straight  sections 
Kl  ~  |  Le I D  for  channel  bend 


(73) 

(74) 


a  cos  0C 


(1.417  —  4.240s  +  3.789s2) 


ds 

dl 


c3dPg  ViPiVii 
S  dl  +  I<p 


s3p,g  cos  6,/g 


(65) 


For  porous  media,  only  liquid-water  transport  in  the  z-direction 
is  considered  and  the  gas-phase  pressure  gradient  in  the  z- 
direction,  dPgldz  is  ignored.  Therefore,  the  related  governing 
equations  can  be  obtained  by  substituting  the  z  coordinate  for  the 
l  coordinate  and  taking  dPg/dz  as  zero. 

For  hydrophilic  porous  media  of  6c  <  90° : 

aC0Sdc(rp 

=  V,P'Vl'z  -  s3  Pig  cos  ez/g  (66) 

KP 


1/2 


s3[-1.417 +  4.240(1  -s)-3.789(l  -s)2]  — 

dz 


It  can  be  determined  from  Eqs.  (73)  and  (74)  that  the  I<p  of  straight 
channel  sections  depends  on  both  D  and  wc/dc,  while  the  I<p  of 
channel  bends  is  related  to  Le/D  as  well  as  D  and  wc/dc.  Therefore, 
it  may  be  predicted  that  compared  with  straight  channel  sections,  a 
smaller  K£  is  formed  for  channel  bends,  and  this  can  lead  to  appar¬ 
ent  differences  of  flooding  within  the  channel. 

In  Eqs.  (64)-(67),  p/V/  /  and  p/V/z  are  variables  that  both  need  to 
be  solved  and  can  be  determined  as  follows: 


dpiVu  dpiVlz 
~dT  ~  dz~  +Ss 


(75) 


dpiVfz  _ 
dz  ~ 


Ss  = 


Mh2o  kc 

^H20^e 


e(l  -S)xy2° 

R7-  (Pg2°-Psat) 

£(M~5)5P'(PgH20-P^) 

Mh20  s 


if  pH 2°  >  Psa[ 
if  pH 2°  <  Psat 


(76) 

(77) 


Eq.  (75)  describes  liquid-phase  mass  conservation  within  the 
oxidant  channel,  in  which  the  liquid-phase  mass  flow  rate  per  unit 
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volume  along  the  streamwise  direction,  dp/V///d/,  is  equal  to  the 
sum  of  the  liquid-phase  mass  flow  rate  per  unit  volume  in  the 
z-direction,  -dp/V/z/dz,  and  the  liquid-phase  source  within  the 
channel,  Ss.  Eq.  (76)  describes  liquid-phase  conservation  within 
porous  media,  where  the  liquid-phase  mass  flow  rate  per  unit 
volume  in  the  z-direction,  dp/V//dz,  is  equal  to  the  liquid-phase 
source  therein,  5S.  Ss  can  be  determined  using  Eq.  (77),  where  the 
first  and  second  terms  represent  the  liquid-phase  source  and  sink, 
respectively,  both  of  which  result  from  phase  changes  between 
water  vapor  and  liquid  water.  In  addition,  kc  and  ke  are  the  con¬ 
densation  and  evaporation  constant,  respectively,  and  is  the 
water-vapor  molar  fraction.  It  can  be  determined  from  Eq.  (77) 
that  when  the  water-vapor  pressure,  P|*2°,  is  larger  than  the  sat¬ 
uration  pressure,  Psau  Ss  represents  the  condensation  rate  from 
water  vapor  to  liquid  water,  while  when  P|b°  is  smaller  than 
Psau  Ss  represents  the  evaporation  rate  from  liquid  water  to  water 
vapor. 


b 


u 


Original  2-channel  flow-field 
including  34  bends 


Original  components  including  the 
CGDL,  CCL,  PEM,  ACL,  and  AGDL 


3.  Results  and  discussion 

To  verify  the  accuracy  and  reliability  of  the  computational  mod¬ 
els  developed  in  this  study,  three  cases  were  analyzed  in  sequence 
to  perform  a  systematic  evaluation  and  examination.  Original  fuel¬ 
cell  models  for  each  case  and  the  simplified  ones  adopted  by  this 
study  are  shown  in  Fig.  3(a)-(c),  respectively,  and  detailed  physi¬ 
cal  parameters  and  operating  conditions  are  tabulated  in  Table  1. 
In  these  cases,  case  1  is  the  fuel-cell  experiment  presented  by 
Ticianelli  et  al.  [36],  and  cases  2  and  3  are  fuel-cell  simulations  pre¬ 
sented  by  Wang  and  Wang  [14,15].  Although  these  studies  give 
detailed  operating  conditions,  they  cannot  provide  all  of  the  phys¬ 
ical  parameters  necessary  for  the  present  computational  models. 
For  unavailable  physical  parameters,  this  study  assumes  reason¬ 
able  values  to  define  them.  The  operating  conditions  and  physical 
parameters  for  these  cases  are  discussed  further  in  the  following 
paragraphs. 


c 


Simplified  components  including 
the  CGDL,  CCL,  and  half  PEM 


Fig.  3.  (a)  Original  fuel-cell  models  and  simplified  ones  for  case  1.  (b)  Original  fuel-cell  models  and  simplified  ones  for  case  2.  (c)  Original  fuel-cell  models  and  simplified 
ones  for  case  3. 
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Table  1 

Physical  parameters  and  operating  conditions  for  each  case. 


Quantity 


Case  1 


Case  2 


Case  3 


Gas  channel  pattern 
Gas  channel  length,  Lf“ 

Gas  channel  depth,  d^as 

Gas  channel  width,  w^as 

Shoulder  width  between  gas  channels,  wfas 

Coolant  channel  pattern 

Coolant  channel  length,  L™olant 

Coolant  channel  depth,  dccoolant 

Coolant  channel  width,  wccoolant 

Shoulder  width  between  coolant  channels, 

yycoolant 

BP  thickness,  dBp 

BP  electrical  resistivity  in  x-y  direction,  RBp,Xy 
BP  electrical  resistivity  in  z  direction,  RBPtZ 
BP  thermal  conductivity,  /<JP 
BP  contact  angle,  Of 
GDL  thickness,  dGD; 

GDL  porosity,  eGDL 

GDL  tortuosity,  tGdl 

GDL  permeability,  I<fiDL 

GDL  electrical  resistivity  in  x-y  direction, 

RGDL,xy 

GDL  electrical  resistivity  in  z  direction,  RGdl,z 

GDL  thermal  conductivity,  kfDL 

GDL  contact  angle,  6fDL 

CL  thickness,  da 

CL  porosity,  sGl 

PEM  thickness,  dpBM 

PEM  porosity,  sPEm 

Contact  electrical  resistance  between  the  rib 
and  GDL,  R^GDL 

Contact  thermal  resistance  between  the  rib 
and  GDL,  Rf /GDL 

Dimensionless  equivalent  length,  Le/D 
Density  of  liquid  water,  pi 
Kinematic  viscosity  of  liquid  water,  v/ 
Thermal  conductivity  of  liquid  water,  Iq 
Surface  tension,  cr 
Gravity,  g 

Molecular  weight  of  N2,  MNl 
Molecular  weight  of  02,  Mo2 
Molecular  weight  of  H20,  Mh2o 
Experimental  constant,  a  for  N2  and  02/H20 
Experimental  constant,  b  for  N2  and  02/H20 
Critical  pressure,  pci  for  N2/02/H20 
Critical  temperature,  Tci  for  N2/02/H20 
Experimental  constant,  A  for  N2/02/H20 

Experimental  constant,  B,  for  N2/02/H20 
Experimental  constant,  a*j  for  N2/02/H20 
Experimental  constant,  al2  for  N2/02/H20 


Straight  pattern 

71.12  mm 

0.762  mm 

0.762  mm 

0.762  mm 

NA 

NA 

NA 

NA 

NA 

2.286  mm 
5  x  10-5  £2  m 
24  x  10“5  £2  m 
20  Wm1  k"1 
80° 

0.254  mm 

0.4 

1.5 

1.76  x  10-11  m2 
0 

5.77  x  10-5  £2  m 
1.7  Wm-1  k"1 
100° 

0.0287  mm 
0.4 

0.23  mm 

0.28 

0 

0 

0 

977.8  kg  m~3 
4.14  x  10“7  m2  s_1 
0.66  Wm-1  k-1 
0.0625  Nm-1 

9.8  ms-2 
28  g  moD1 
32  g  mol-1 
18  g  mol-1 

2.745  x  10-4/3.640x  10“4 
1.823/2.334 
33.5/49.7/221.2  atm 
126.2/154.4/647.3 1< 

1.405  x  10-6/1.783  x  10-6/1.867  x  10“6  kgs"1 

m-1 1<-°-5 

111.5/156/708  K 

3.298677/3.212936/3.386842 

1.40824  x  10_3/1. 127486  x  10-3/3.474982 

x  10-3  K1 


Serpentine  pattern 

Serpentine  pattern 

1299  mm 

431  mm 

1.0  mm 

1.0  mm 

1.0  mm 

1.0  mm 

1.0  mm 

1.0  mm 

NA 

Serpentine  pattern 

NA 

431  mm 

NA 

1.0  mm 

NA 

1.0  mm 

NA 

1.0  mm 

3.0  mm 

3.0  mm 

5  x  10-5  £2  m 

5  x  10-5  £2  m 

24  x  10-5  £2  m 

24  x  10-5  £2  m 

20  Wm-1  k"1 

20  Wm-1  k"1 

80° 

80° 

0.3  mm 

0.2  mm 

0.6 

0.6 

1.5 

1.5 

10“12  m2 

10“12  m2 

0 

0 

5.77  x  10-5  £2  m 

5.77  x  10-5  £2  m 

1.7  Wm-1  k"1 

1.7  Wm-1  k-1 

100° 

100° 

0.01  mm 

0.01  mm 

0.4 

0.4 

0.051  mm 

0.045  mm 

0.28 

0.28 

0 

0 

0 

0 

60 

60 

Experimental  constant,  a '3  for  N2/02/H20 
Experimental  constant,  a\  for  N2/02/H20 
Experimental  constant,  a'5  for  N2/02/H20 


-3.963222  x  10-6/-5.75615  x 
10_7/-6.354696  x  10-6  K-2 
5.641515  x  10-9/1.313877x 
10-9/6.968581  x  1 0  9  K  3 
-2.444855  x  10-12/-8.768554  x 
10_13/-2. 506588  x  10“12  I<“4 


Universal  ideal  gas  constant,  R 

8.3145  kjkmol"1  K_1 

Faraday  constant,  F 

96.485C  mol-1 

Condensation  rate  constant,  kc 

100  s-1 

Evaporation  rate  constant,  ke 

1  atm-1  s_1 

Reference  pressure,  Prgef 

101,325  Pa 

Activation  energy,  Ec 

66  kj  moD1 

Reference  temperature,  Tref 

298.15  K 

Transfer  coefficient,  a 

0.5 

Cathode  reference  exchange  current  density, 

?ef 

0.1 5  Am-2 

0.20  Am-2 

0.20  Am-2 

lo 

Fuel/oxidant 

H2/air 

H2/air 

H2/air 

Anode/cathode  inlet  flow  rate,  N™^a/Nl^c 

2.8  x  104/3.0  x  104 

1?g<  103A  m-2: 

i7<103Am-2: 

Am-2  equivalent 

2.0  x  103/2.0  x  103 

2.0  x  103/2.0  x  103 

Am-2  equivalent 

Am-2  equivalent 

izavg  >  103A  m-2: 

izavg  >  103A  m-2: 

2.0/2.0  stoichiometric 

2.0/2.0  stoichiometric 

Anode/cathode  inlet  temperature,  Vf/Tf1 

50/50  °C; 

80/80  °C 

80/80  °C 

80/80  °C 
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Table  1  ( Continued ) 


Quantity 

Case  1 

Case  2 

Case  3 

Anode/cathode  inlet  relative  humidity, 

100%/100% 

100%/0% 

50%/50% 

RHa  /RHC 

in'  in 

Anode/cathode  outlet  pressure,  PfXt,a/ptg°out,c 

l.O/l.Oatm; 

3.0/5.0atm 

2.0/2.0atm 

2.0/2.0atm 

Fuel  cell  temperature,  Tceu 

50°C; 

80°C 

80°C; 

80°C 

72-82  °C  with  the 
oxidant  flow 

First,  in  case  1,  experimental  polarization  curves  were  mea¬ 
sured  with  fixed  reactant  flow  rates,  where  the  anode  and  cathode 
inlet  flow  rates,  N|^a  and  N^c,  are  equivalent  to  those  generat¬ 
ing  current  densities  of  2.8  x  104  and  3.0  x  104  Am-2,  respectively, 
and  the  oxidant  and  fuel  flows  are  parallel.  Under  fixed  flow 
rates,  three  different  operating  conditions  were  used,  where  the 
anode/cathode  inlet  temperature  and  the  outlet  pressure  are 
50/50°C  and  l.O/l.Oatm,  50/50°C  and  3.0/5.0atm,  and  80/80°C 
and  3.0/5.0atm,  respectively.  In  this  study,  these  three  operating 
conditions  are  represented  as  case  1(a),  case  1(b)  and  case  1(c), 
respectively,  for  convenience.  In  the  experiments  by  Ticianelli  et  al. 
[36],  how  the  fuel  cell  is  oriented  relative  to  gravity  is  unknown.  In 
order  to  obtain  6[/g  and  0z/g ,  which  are  necessary  for  the  compu¬ 
tational  models,  it  is  assumed  that  fluid  flows  within  the  oxidant 
channel  cannot  be  affected  by  gravity,  while  fluid  flows  within  the 
CGDL  are  parallel  to  gravity.  Therefore,  &i/g  and  6z/g  are  specified  as 
90°  and  0°,  respectively.  In  addition,  for  predictions  of  polarization 
curves  of  case  1,  this  study  chooses  13  operating  points,  including 
1 0, 20, 50,  and  1 00-1000  mA  cm-2  with  an  interval  of  1 00  mA  cm-2, 
to  simulate  cell  voltages  corresponding  to  specific  current  densi¬ 
ties.  From  the  simulation  of  case  1  performed  in  this  study,  the 
predicted  and  experimental  polarization  curves  are  compared  to 
evaluate  whether  computational  models  can  yield  qualitative  and 
quantitative  agreement  with  experiments  under  different  operat¬ 
ing  conditions. 

In  case  2,  a  3D  fuel-cell  model,  which  includes  the  reactant  flow 
fields  of  the  anode  and  cathode,  was  built  by  Wang  and  Wang  [14], 
and  in  case  3,  a  3D  fuel-cell  model  that  includes  the  reactant  flow 
fields  of  the  anode  and  cathode  and  the  coolant  flow  fields  was 
built  by  Wang  and  Wang  [15].  For  cases  2  and  3,  the  active  areas 
are  50  and  200  cm2,  respectively,  and  the  flow  fields  are  a  two- 
channel  serpentine  pattern  with  a  channel  length  of  1299  mm  and 
a  24-channel  serpentine  pattern  with  a  channel  length  of  431  mm, 
respectively.  Case  2  is  simulated  numerically  for  cell  voltages  of 
0.65  V  under  specific  operating  conditions,  and  case  3  for  cur¬ 
rent  densities  of  1 A  cm-2  under  two  different  cell  temperatures, 
which  are  a  constant  80  °C  and  an  increasing  72-82  °C  with  oxi¬ 
dant  flows,  respectively.  Therefore,  the  predictions  of  Wang  and 
Wang  [14,15]  include  only  a  single  point  of  current  density  vs.  cell 
voltage  (IV);  therefore,  no  polarization  curves  are  available  for  com¬ 
parison.  In  addition,  the  reactant  flow  rates  of  both  cases  are  given 
with  stoichiometric  feeding,  where  N[^a  and  Nl^c  are  equivalent 
to  stoichiometric  flow  rates  of  2.0  and  2.0,  respectively.  However, 
in  this  study,  two-step  flow  rates  are  used  to  simultaneously  sat¬ 
isfy  reactant  feeding  for  the  simulated  points  of  Wang  and  Wang 
[14,15]  and  stability  and  convergence  in  numerical  computations. 
That  is,  when  the  average  current  density,  izvg,  is  inferior  or  equal  to 
103  Am-2,  N[^a  and  Nl^c  are  both  equivalent  to  flow  rates  generat¬ 
ing  current  densities  of  2.0  x  103  Am-2,  while  when  izvg  is  greater 
than  1 03  A  nrr 2 ,  N[^  a  and  N ™ ’ c  are  equivalent  to  stoichiometric  flow 
rates  of  2.0  and  2.0,  respectively. 

In  the  research  of  Wang  and  Wang  [14,15],  it  is  also  unknown 
how  the  fuel  cell  is  oriented  relative  to  gravity.  Therefore,  6yg  and 
0Zjg  are  similarly  specified  as  90°  and  0°,  respectively,  to  give  the 
physical  parameters  necessary  for  cases  2  and  3.  For  predictions 


of  the  polarization  curves  of  cases  2  and  3,  this  study  chooses  18 
operating  points,  including  1 0, 20, 50,  and  1 00-1500  mA  cm-2  with 
an  interval  of  100  mA cm-2,  to  simulate  cell  voltages  correspond¬ 
ing  to  specific  current  densities.  In  case  3,  simulations  under  a 
constant  80  °C  or  an  increasing  72-82  °C  with  oxidant  flows  are  rep¬ 
resented  as  case  3(a)  and  case  3(b),  respectively.  The  simulations 
of  cases  2  and  3  performed  in  this  study  are  focused  on  evaluating 
whether  computational  models  can  predict  overall  performance 
characteristics,  such  as  polarization  curves,  reactant  flow  rates  and 
pressure  drops,  in  an  efficient  manner  and  capture  the  IV  points 
predicted  by  Wang  and  Wang  [  1 4,1 5  ].  Therefore,  the  computational 
resources  and  the  time  required  between  this  study  and  Wang  and 
Wang  [  1 4,1 5  ]  will  also  be  compared.  The  results  and  discussion  are 
divided  into  five  sections,  which  are  comparisons  between  exper¬ 
iments  and  predictions  for  case  1,  predicted  open-circuit  voltages 
and  polarization  losses  of  case  1,  the  predicted  overall  performance 
of  case  2,  the  predicted  overall  performance  of  case  3,  and  compar¬ 
isons  of  modeling  ability  between  different  computational  models. 

3.1.  Comparisons  between  experiments  and  predictions  for  case  I 

Figs.  4-6  show  the  comparisons  between  experiments  and  pre¬ 
dictions  for  case  1  under  different  operating  conditions.  It  can 
be  determined  by  observing  the  overall  polarization  curves  from 
Figs.  4-6  that  numerical  predictions  indeed  can  capture  exper¬ 
imental  data  varying  with  different  operating  conditions,  both 
qualitatively  and  quantitatively.  First,  for  case  1(a),  with  oper¬ 
ating  conditions  of  50/50 °C  and  l.O/l.Oatm,  the  experimental 
and  predicted  /V  values  vary  from  1.000  V  at  0  mA  cm-2  and 
1.201  V  at  0  mA  cm-2  to  0.050  V  at  1000  mA  cm-2  and  0.331  V  at 
1000  mA  cm-2,  respectively.  Then,  for  case  1(b),  with  operating 
conditions  of  50/50  °C  and  3.0/5.0  atm,  experimental  and  predicted 
IV  values  vary  from  1 .020  V  at  0  mA  cm-2  and  1 .228  V  at  0  mA  cm-2 
to  0.350  V  at  1000  mA cm-2  and  0.413  V  at  1000  mA  cm-2,  respec¬ 
tively.  Finally,  for  case  1(c)  with  operating  conditions  of  80/80  °C 
and  3.0/5.0atm,  experimental  and  predicted  IV  values  vary  from 
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Fig.  4.  Comparisons  between  experiments  and  predictions  for  case  1(a). 
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Fig.  5.  Comparisons  between  experiments  and  predictions  for  case  1(b). 
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Fig.  7.  Predicted  open-circuit  voltages  and  polarization  losses  for  case  1(a). 
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Fig.  6.  Comparisons  between  experiments  and  predictions  for  case  1(c). 

1.020  V  at  0  mA  cm-2  and  1.204  V  at  0  mA  cm-2  to  0.520  V  at 
1000mA cm-2  and  0.535 Vat  1000mA cm-2,  respectively. 

Further  comparisons  of  the  experimental  and  predicted  IV 
values  for  case  1(a)  reveal  that  the  predicted  cell  voltages  under¬ 
estimate  the  experimental  data  for  low  current  densities  between 
0  and  200  mA  cm-2  and  intermediate  current  densities  between 
200  and  600mA cm-2,  while  they  overestimate  the  experimental 
data  for  high  current  densities  between  600  and  1000  mA cm-2. 
For  case  1(b),  the  predicted  cell  voltages  fit  the  experimental  data 
for  low  current  densities  between  0  and  200  mA  cm-2  and  inter¬ 
mediate  current  densities  between  200  and  600  mA  cm-2,  while 
they  overestimate  the  experimental  data  for  high  current  densi¬ 
ties  between  600  and  1000  mA cm-2.  For  case  1(c),  the  predicted 
cell  voltages  overestimate  the  experimental  data  for  low  current 
densities  between  0  and  200  mA cm-2,  while  they  fit  the  exper¬ 
imental  data  for  intermediate  current  densities  between  200  and 
600  mA  cm-2  and  high  current  densities  between  600  and  1 000  mA 
cm-2. 

The  reasons  for  the  differences  between  the  experiments  and 
the  predictions  may  include  the  fact  that  computational  models 
cannot  describe  physical  phenomena  precisely  and  the  fact  that 
the  operating  conditions  cannot  be  controlled  accurately;  hence, 
the  computational  models  use  incorrect  parameters.  Although 
Ticianelli  et  al.  [36]  provided  only  one  set  of  experimental  param¬ 
eters  and  no  detailed  flow  rates  and  cell  temperatures  under 
different  current  densities,  the  computational  models  developed 
by  this  study  can  still  predict  the  changing  trends  of  polarization 
curves  corresponding  to  experiments  under  different  operating 
conditions.  Therefore,  from  the  comparisons  between  predictions 
and  experiments  for  case  1,  the  accuracy  and  reliability  of  the 


Fig.  8.  Predicted  open-circuit  voltages  and  polarization  losses  for  case  1(b). 


computational  models  is  demonstrated.  In  order  to  investigate  the 
effects  of  open-circuit  voltages  and  polarization  losses  on  cell  volt¬ 
ages,  as  shown  in  Eq.  (59),  the  various  components  forming  cell 
voltages  are  described  and  analyzed. 

3.2.  Predicted  open-circuit  voltages  and  polarization  losses  for 
case  I 

Figs.  7-9  show  predicted  open-circuit  voltages  and  polarization 
losses  under  different  operating  conditions  for  case  1,  where  the 
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Fig.  9.  Predicted  open-circuit  voltages  and  polarization  losses  for  case  1(c). 
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x-axis  represents  the  current  density,  and  the  y-axis  represents 
the  various  components  that  constitute  the  cell  voltages.  These 
components  contain  average  open-circuit  voltages,  V£c,  average 
overpotentials,  r)Cl  average  ionic-phase  ohmic  losses,  Sty,  and  aver¬ 
age  electronic-phase  ohmic  losses,  8(j)s.  It  can  be  determined  from 
Eq.  (59)  that  cell  voltages  ty  -  ty  can  be  expressed  as  v£c  -rjc  - 
Sty  -  Sty,  and  thus  average  cell  voltages  ty  -  ty  can  be  expressed 
as  V&  -fjc-  84>s  -  Sty.  According  to  ty  -  ty  =  V£c  -  rjc  -  8(j)s  - 
8(pf ,  variations  of  average  cell  voltages  with  increasing  current  den¬ 
sities  under  different  operating  conditions  can  be  explained  by 
analyzing  each  component  of  the  cell  voltage. 

First,  it  can  be  determined  from  observing  Vfc  in  respective 
figures  that  when  current  densities  increase,  v£c  varies  little.  Com¬ 
paring  v£  under  the  same  current  densities,  it  is  found  that 
V0CC  increases  with  increasing  anode  and  cathode  outlet  pres¬ 
sures.  That  is,  v£  increases  from  1.201V  at  10mA  cm-2-l. 196  V 
at  1000  mA  cm-2  in  Fig.  7  to  1.229  V  at  10  mA  cm-2 -1.224  V  at 
1 000  mA  cm-2  in  Fig.  8.  However,  when  cell  temperatures  increase, 
V0CC  decreases  instead,  where  v£c  decreases  from  the  values  shown 
in  Fig.  8  to  1 .204  V  at  1 0  mA  cm-2-l  .200  V  at  1 000  mA  cm-2  in  Fig.  9. 
Then,  by  observing  rjc,  it  is  determined  that  this  value  increases 
with  increasing  current  densities.  By  comparing  rjc  under  the  same 
current  densities,  it  is  determined  that  rjc  decreases  with  increasing 
anode  and  cathode  outlet  pressures  and  cell  temperatures.  That  is, 
r)c  decreases  from  0.295  V  at  1 0  mA  cm-2  -0.560  V  at  1 000  mA  cm-2 
in  Fig.  7  to  0.248  V  at  1 0  mA  cm-2 -0.51 2  V  at  1 000  mA  cm-2  in  Fig.  8, 
and  then  to  0.147  V  at  10  mA cm-2 -0.436  V  at  1000  mA  cm-2  in 
Fig.  9.  __ 

Additionally,  it  is  observed  that  Sty  increases  with  increasing 
current  densities.  When  comparing  Sty  under  the  same  cur¬ 
rent  densities,  it  can  be  determined  that  Sty  decreases  with 
increasing  anode  and  cathode  outlet  pressures  and  cell  tempera¬ 
tures.  That  is,  Sty  decreases  from  0.003  V  at  10  mA  cm-2 -0.266  V 
at  1000  mA  cm-2  in  Fig.  7  to  0.003  V  at  10  mA  cm-2 -0.259  V  at 
1 000  mA  cm-2  in  Fig.  8,  and  then  to  0.002  V  at  1 0  mA  cm-2  -0.1 89  V 
at  1000  mA  cm-2  in  Fig.  9.  Finally,  observations  reveal  that  8(j)s 
increases  with  increasing  current  densities.  By  comparing  S(j)s 
under  the  same  current  densities,  it  can  be  determined  that  84>s 
is  consistent;  the  values  of  8(j)s  in  Figs.  7-9  all  vary  from  4  x  10-4  V 
at  10  mA  cm-2  to  0.039  V  at  1000  mA  cm-2. 

Based  on  the  above  description  for  case  1,  overpotentials  are 
the  primary  polarization  losses,  ionic-phase  ohmic  losses  are  the 
secondary,  and  electronic-phase  ohmic  losses  are  the  minimum. 
However,  when  operating  conditions  change,  the  variability  in 
open-circuit  voltages  and  polarization  losses  also  changes.  To 
explain  these  phenomena,  the  properties  governing  open-circuit 
voltages  and  polarization  losses  need  to  be  identified  first  so  that 
the  interactions  between  various  components  of  cell  voltages  and 
governing  properties  can  be  further  understood.  The  distributions 
of  governing  properties  will  be  presented  in  the  second  part  of 
this  series  of  studies  to  analyze  and  discuss  the  effects  of  these 
properties  on  various  components  of  cell  voltages  in  depth. 

3.3.  Modeling  results  of  overall  performance  characteristics  for 
case  2 

Figs.  10-12  show  the  comparisons  of  predicted  IV  values  and 
output  powers  between  different  computational  models,  predicted 
oxidant  flow  rates  and  pressure  drops,  and  predicted  open-circuit 
voltages  and  polarization  losses,  respectively,  for  case  2.  Fig.  10 
shows  that  the  IV  values  predicted  by  the  computational  mod¬ 
els  of  this  study  vary  from  1.191V  at  0  mA  cm-2  to  0.536  V  at 
1500  mA  cm-2  and  that  the  predicted  powers  vary  from  0W  at 
0  mA  cm-2  to  40.2  W  at  1 500  mA  cm-2 .  For  a  cell  voltage  of  0.64  V, 
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Fig.  10.  Comparisons  of  predicted  IV  values  and  output  powers  between  different 
computational  models  for  case  2. 
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Fig.  11.  Predicted  oxidant  flow  rates  and  pressure  drops  for  case  2. 


the  predicted  current  density  and  the  power  are  900  mAcm-2  and 
28.8  W,  respectively.  For  simulations  at  a  cell  voltage  of  0.65  V, 
performed  by  Wang  and  Wang  [14],  considering  and  ignoring 
convection  effects  within  GDLs  can  result  in  different  current  den¬ 
sities  and  output  powers,  which  are  880 mAcm-2  and  28.6  W,  and 
910 mAcm-2  and  29.6 W,  respectively.  Therefore,  compared  with 
the  predictions  of  Wang  and  Wang  [14],  the  computational  models 
of  this  study  can  not  only  accurately  reproduce  the  current  density 
and  output  power  at  specific  cell  voltage,  but  can  also  quickly  give 
predictions  at  other  cell  voltages  to  produce  polarization  curves  and 
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Fig.  12.  Predicted  open-circuit  voltages  and  polarization  losses  for  case  2. 
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Fig.  13.  Comparisons  of  predicted  IV  values  and  output  powers  between  different 
computational  models  for  case  3. 


corresponding  powers,  which  are  not  available  from  the  original 
research  [14]. 

Fig.  1 1  shows  that  the  predicted  oxidant  flow  rates  vary  from 
Oslpm  at  0  mA  cm-2  to  2.633  slpm  at  1500  mA  cm-2  and  that  the 
predicted  pressure  drops  vary  from  0  Pa  at  0  mA  cm-2  to  40,264  Pa 
at  1500mA cm-2.  At  a  current  density  of  900mA cm-2,  the  pre¬ 
dicted  oxidant  flow  rates  and  pressure  drops  are  1.580  slpm  and 
24,630  Pa,  respectively.  Finally,  it  can  be  determined  from  Fig.  12 
that  when  current  densities  increase,  V£c  varies  little  with  the 
value  ranging  between  1.191V _at  10  mA  cm-2  and  1.183  V  at 
1 500  mA  cm-2 ,  while  rjc ,  Sfy  and  all  increase.  By  comparing  var¬ 
ious  polarization  losses,  it  is  found  that  that  rjc  is  the  primary  loss, 
with  a  value  ranging  between  0.154  V  at  10  mA  cm-2  and  0.466  V 
at  1500mA cm-2;  Sfy  is  the  secondary  loss,  with  a  value  ranging 
between  0.001  V  at  1 0  mA  cm-2  and  0.1 1 0  V  at  1 500  mA  cm-2 ;  and 
8(j)s  is  the  minimum  loss,  with  a  value  ranging  between  5  x  10-4  V 
at  10  mA cm-2  and  0.070  V  at  1500  mA cm-2. 

3.4.  Modeling  results  of  overall  performance  characteristics  for 
case  3 

Figs.  13-16  show  the  comparisons  of  predicted  IV  values  and 
output  powers  between  different  computational  models,  predicted 
oxidant  flow  rates  and  pressure  drops,  and  predicted  open-circuit 
voltages  and  polarization  losses,  respectively,  for  case  3.  Fig.  13 
shows  that  for  case  3(a),  with  a  cell  temperature  of  a  constant 
80  °C,  the  IV  values  predicted  by  the  computational  models  in  this 
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Fig.  15.  Predicted  open-circuit  voltages  and  polarization  losses  for  case  3(a). 


study  vary  from  1.190  V  at  0  mA  cm-2  to  0.481  V  at  1500  mA  cm-2. 
The  predicted  powers  vary  from  0W  at  0  mA  cm-2  to  144.3W  at 
1500  mA  cm-2,  where  at  a  current  density  of  1000  mA  cm-2,  the 
predicted  cell  voltage  and  power  are  0.590  V  and  118.0W,  respec¬ 
tively.  For  case  3(b),  with  a  cell  temperature  increasing  with  oxidant 
flows  from  72  to  82  °C,  the  predicted  /V  values  vary  from  1 .1 93  V  at 
0  mA  cm-2  to  0.542  V  at  1 500  mA  cm-2.  The  predicted  powers  vary 
from  0  W  at  0  mA  cm-2  to  1 62.6  W  at  1 500  mA  cm-2,  where  at  a  cur¬ 
rent  density  of  1 000  mA  cm-2,  the  predicted  cell  voltage  and  power 
are  0.624  V  and  124.8  W,  respectively.  In  simulations  for  a  current 
density  of  1000  mA  cm-2  performed  by  Wang  and  Wang  [15],  case 
3(a)  and  case  3(b)  can  result  in  different  cell  voltages  and  output 
powers,  which  are  0.590  V  and  118.0W  and  0.616  V  and  123.2  W, 
respectively.  For  case  3,  the  computational  models  in  this  study 
can  not  only  reproduce  cell  voltage  and  output  power  at  a  specific 
current  density,  but  can  also  predict  polarization  curves  and  corre¬ 
sponding  powers  that  cannot  be  achieved  by  the  model  developed 
in  the  original  research  [15]. 

Fig.  14  shows  that  for  case  3(a)  and  case  3(b),  the  predicted 
oxidant  flow  rates  vary  from  Oslpm  at  0mA cm-2  to  11.884 slpm 
at  1 500  mA  cm-2  and  from  0  slpm  at  0  mA  cm-2  to  1 1 .882  slpm  at 
1500  mA cm-2,  respectively.  Furthermore,  the  predicted  pressure 
drops  vary  from  0  Pa  at  0  mA  cm-2  to  2280  Pa  at  1 500  mA  cm-2  and 
from  OPa  at  0  mA cm-2  to  2583  Pa  at  1500  mA cm-2,  respectively. 
In  addition,  at  a  current  density  of  1000  mA  cm-2,  the  predicted 
oxidant  flow  rate  and  pressure  drop  are  7.926  slpm  and  1528  Pa, 
respectively,  for  case  3(a)  and  7.925  slpm  and  1 714  Pa,  respectively, 
for  case  3(b).  Finally,  it  can  be  determined  from  Figs.  15  and  16 


Current  density  (mA  cm'2) 

Fig.  14.  Predicted  oxidant  flow  rates  and  pressure  drops  for  case  3. 


Current  density  (mA  cm'2) 


Fig.  16.  Predicted  open-circuit  voltages  and  polarization  losses  for  case  3(b). 
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Table  2 

Comparisons  of  computational  resources  and  modeling  times  between  different  computational  models  for  case  2. 
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Case  2 

Computational  resources 

Modeling  time 

The  model  used  by  Wang  and  Wang  [14] 

9  nodes  of  1.4  GHz  AMD  Athelon  Thunderbird  CPU  and 

512  MB  DDR  SDRAM 

Nearly  5  h  for  only  one  IV  point 

The  model  developed  in  this  study 

2.39  GHz  Dual  Core  AMD  Opteron  CPU  and  3.5  GB  DDR  SDRAM 

Nearly  4  h  for  a  polarization  curve  including  18  IV  points 

Table  3 

Comparisons  of  computational  resources  and  modeling  times  between  different  computational  models  for  case  3. 

Case  3 

Computational  resources 

Modeling  time 

The  model  used  by  Wang  and  Wang  [15] 

32  nodes  of  2.8  GHz  Intel  Pentium  4  CPU  and  1.0  GB  DDR  SDRAM 

Nearly  20  h  for  only  one  IV  point 

The  model  developed  in  this  study 

2.39  GHz  Dual  Core  AMD  Opteron  CPU  and  3.5  GB  DDR 
SDRAM 

Nearly  1  and  2  h  for  a  polarization  curve  including  18 

IV  points  of  case  3(a)  and  case  3(b),  respectively 

that  when  current  densities  increase,  the  v£c  of  both  cases  3(a) 
and  3(b)  vary  little,  the  values  of  which  range  between  1.190  V 
at  10  mA  cm-2  and  1.184  V  at  1500  mA  cm-2  and  between  1.193  V 
at  10  mA cm-2  and  1.186  V  at  1500  mA cm-2,  respectively.  How¬ 
ever,  with  increasing  current  densities,  rjc ,  Sty  and  8(j)s  all  increase 
for  cases  3(a)  and  3(b).  Comparing  various  polarization  losses  in 
Figs.  15  and  16,  it  is  found  that  rjc  increases  from  0.158  V  at 
1 0  mA cm-2 -0.467  V  at  1500  mA  cm-2  for  case  3(a)_to  0.168  V  at 
1 0  mA cm-2 -0.479  V  at  1500  mA  cm-2  for  case  3(b).  8ty  decreases 
from  0.002  Vat  10  mA  cm-2 -0.1 68  Vat  1500  mA  cm-2  for  case  3(a) 
to  0.001  V  at  10  mA  cm-2 -0.098  V  at  1500  mA  cm-2  for  case  3(b). 
8(j)s  is  the  same  for  case  3(a)  and  case  3(b)  with  a  value  ranging 
between  5x1 0-4  V  at  1 0  mA  cm-2  and  0.068  V  at  1 500  mA  cm-2. 

Based  on  the  above  descriptions,  rjc  is  the  primary  polarization 
loss,  Sty  is  the  secondary,  and  is  the  minimum.  By  observing 
further  that  performance  differences  between  different  cell  tem¬ 
peratures  are  dominated  by  polarization  loss,  it  can  be  determined 
that  under  the  same  current  density,  the  overpotentials  of  case  3(a) 
are  slightly  higher  than  those  in  case  3(b),  while  the  ionic-phase 
ohmic  losses  of  case  3(a)  are  much  lower  than  those  of  case  3(b). 
This  leads  to  the  finding  that  the  total  polarization  losses  of  case 
3(b)  are  smaller  than  for  case  3(a)  because  of  smaller  ionic-phase 
ohmic  losses  caused  by  lower  cell  temperatures. 

3.5.  The  comparisons  of  modeling  ability  between  different 
computational  models  for  cases  2  and  3 

In  order  to  compare  the  modeling  ability  of  the  computa¬ 
tional  models  in  this  study  to  that  of  Wang  and  Wang  [14,15], 
Tables  2  and  3  give  computational  resources  and  modeling  times 
for  cases  2  and  3,  respectively.  First,  Table  2  tabulates  that  in  case 
2,  Wang  and  Wang  [14]  takes  nearly  5  h  to  implement  a  simulation 
for  only  one  IV  point  on  9  nodes  of  a  1 .4  GHz  AMD  Athelon  Thunder- 
bird  CPU  with  512  MB  DDR  SDRAM.  This  study  takes  almost  4  h  to 
implement  18  simulations  for  the  polarization  curve,  including  18 
IV  points  on  the  HP  xw9300  Workstation  of  one  node  of  a  2.39  GHz 
Dual  Core  AMD  Opteron  CPU  with  3.5  GB  DDR  SDRAM.  Second, 
Table  3  tabulates  that  in  case  3,  Wang  and  Wang  [15]  takes  nearly 
20  h  to  implement  a  simulation  for  only  one  IV point  on  32  nodes  of 
a  2.8  GHz  Intel  Pentium  4  CPU  with  1 .0  GB  DDR  SDRAM.  This  study 
takes  almost  1  and  2h,  respectively,  to  implement  18  simulations 
for  the  polarization  curve,  including  18  IV  points  of  case  3(a)  and 
case  3(b)  on  the  same  HP  xw9300  Workstation. 

Because  this  study  adopts  the  models  covered  by  one  single 
oxidant  channel  as  computational  domains,  case  3,  with  shorter 
total  channel  lengths,  instead  requires  fewer  computational  grids 
relative  to  case  2,  and  this  leads  to  case  3  requiring  less  model¬ 
ing  time  than  case  2.  If  comparing  the  computational  models  of 
this  study  to  that  of  Wang  and  Wang  [14,15],  then  this  study  can 
implement  extensive  modeling  works  in  a  more  efficient  manner 
by  using  fewer  computational  resources.  This  is  because  in  this 


study,  a  cathode  physical  model  that  only  includes  the  regions  cov¬ 
ered  by  one  oxidant  channel  is  considered.  Therefore,  compared 
with  CFD  modeling  with  fuel-cell  models  including  the  regions  cov¬ 
ered  by  all  channels,  great  simplification  of  computational  models 
with  this  study  can  improve  computational  efficiency,  especially 
for  large-scale  fuel-cell  modeling.  In  addition,  this  model,  which 
consists  of  one  dimension  along  the  direction  of  oxidant  flow  and 
one  dimension  along  the  z-direction,  forms  the  main  framework  for 
developing  the  mathematical  formulations.  This  leads  to  the  fact 
that  the  number  of  governing  equations  and  variables  that  need 
to  be  solved  in  this  study  are  much  less  than  those  that  are  solved 
in  three  dimensions;  hence,  computational  efficiency  is  promoted 
again.  Therefore,  combinations  of  simplified  models  and  reduced 
dimensions  are  the  key  to  the  computational  efficiency  of  this 
study,  making  it  superior  to  traditional  CFD  modeling.  This  explains 
why  extensive  modeling  work  can  be  completed  with  fewer  com¬ 
putational  resources  and  in  less  time. 

For  original  fuel-cell  modeling  that  requires  several  hours  to 
several  tens  of  hours,  this  study  presents  refined  and  powerful 
computational  models  as  a  solution  that  can  help  researchers  deal 
with  this  troublesome  issue.  This  study  can  not  only  evaluate  over¬ 
all  performance  characteristics,  such  as  polarization  curves,  output 
powers,  reactant  flow  rates  and  pressure  drops,  but  can  also  provide 
related  property  distributions  to  elucidate  and  explain  microscopic 
phenomena  and  physical  mechanisms  within  fuel  cells.  For  the 
effect  of  property  distributions  on  the  cell  performance  of  the  afore¬ 
mentioned  case  1,  further  analyses  and  discussions  will  be  given  in 
the  second  part  of  this  series  of  studies. 

4.  Conclusions 

This  study  presents  the  model  development  and  validation  of 
an  efficient  method  of  numerical  predictions  for  the  performance 
characteristics  of  fuel  cells.  In  order  to  verify  the  accuracy  and  reli¬ 
ability  of  the  computational  models  in  this  study,  three  cases  were 
analyzed,  and  some  key  conclusions  are  summarized  as  follows: 

1.  From  the  comparisons  between  experimental  and  predicted 
polarization  curves  for  case  1,  the  computational  models  in  this 
study  can  predict  cell  performances  under  different  operating 
conditions,  both  qualitatively  and  quantitatively.  For  cases  2  and 
3,  the  computational  models  of  this  study  can  not  only  accurately 
reproduce  the  IV  value  and  output  power  at  a  specific  operating 
point,  but  can  also  quickly  give  predictions  at  other  operating 
points  to  produce  polarization  curves  and  the  corresponding 
powers. 

2.  In  addition  to  the  polarization  curves  and  output  powers,  the 
computational  models  in  this  study  can  also  predict  overall 
performance  characteristics,  such  as  reactant  flow  rates  and 
pressure  drops,  and  various  components  of  cell  voltages.  These 
modeling  results  can  not  only  evaluate  the  effects  of  flow 
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designs,  operating  conditions  or  material  parameters  on  cell 
performance,  but  can  also  provide  the  reactant  flow  rates  and 
pressure  drops  necessary  for  fuel  cells  so  that  researchers  can 
choose  appropriate  system  components,  such  as  air  blowers. 
Computation  of  open-circuit  voltages  and  various  polarization 
losses  help  researchers  identify  major  factors  dominating  cell 
performance,  which  can  be  further  improved  and  promoted 
according  to  in-depth  analyses  of  various  components  of  cell 
voltages. 

3.  Compared  with  known  simulations  of  large-scale  fuel  cells,  the 
computational  models  of  this  study  can  complete  extensive 
modeling  works  in  a  more  efficient  manner  by  using  fewer  com¬ 
putational  resources.  Therefore,  for  realistic  engineering  designs, 
this  study  presents  a  better  solution,  which  can  not  only  allow 
researchers  to  perform  simulations  of  fuel  cells  at  various  scales 
with  fewer  computational  resources  but  can  also  substantially 
decrease  costs  and  the  time  necessary  for  fuel-cell  experiments. 

For  theoretically  fundamental  investigations,  the  computational 
models  of  this  study  can  also  provide  data  or  pictures  of  property 
distributions.  Therefore,  to  further  investigate  the  effects  of  differ¬ 
ent  operating  conditions  on  cell  performance,  the  distributions  of 
various  components  of  cell  voltages  and  related  properties  govern¬ 
ing  these  components  for  case  1  will  be  given  in  the  second  part  of 
this  series  of  studies.  To  explain  the  interactions  between  various 
properties  and  physical  mechanisms  governing  these  phenomena, 
a  systematical  method  is  used  in  the  second  part  to  perform  detailed 
and  in-depth  analyses. 
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